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We use a singular perturbation method to study the interface dynamics of a non-conserved order param- 
eter (NCOP) system, of the reaction-diffusion type, for the case where an external bias field or convection 
is present. We find that this method, developed by Kawasaki, Yalabik and Gunton for the time-dependant 
Ginzburg-Landau equation and used successfully on other NCOP systems, breaks down for our system when 
the strength of bias/convection gets large enough. 
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I. INTRODUCTION 

The study of non-equilibrium model systems has led to 
much progress in our understanding of self-organisation. 
This organisation emerges when the system is quenched 
to a region of its parameter space where competing mech- 
anisms lead it to a new configuration, which may or may 
not be homogeneous in space or steady in time. The 
transition process often causes interfaces to appear be- 
tween the allowed homogeneous states of the new config- 
uration. These interfaces are diffuse on a molecular scale 
but appear as discontinuities on the mesoscopic length 
scale of the growing domains. The two disparate length 
scales make modelling of such systems and extraction of 
analytical results challenging. Here we will consider a 
model where mechanical transport competes with diffu- 
sive transport and dissipation. 

An analytic method which approximately solves the 
nonlinear partial differential equations (PDEs) modelling 
some pattern forming systems, is a Singular Perturba- 
tion Method (SPM) developed by Kawasaki, Yalabik 
and Gunton (KYG) for the time-dependent Ginzburg- 
Landau equation [1]. It was subsequently applied to 
many other systems where the order parameter is non- 
conserved (NCOP) [1-3], conserved (COP) [4], coupled 
with long range repulsive interactions [5], etc, and which 
exhibit strong nonlinear behaviour at late times. How- 
ever, the KYG method has not yet been applied, to 
our knowledge, to a system with a nonlinear convective 
derivative of the form [v ■ V)'0", where tp is the relevant 
order parameter field in the system. This type of nonlin- 
earity occurs in hydrodynamics and also occurs when an 
external bias field E is present [6]. Then E plays the role 
of V. Hence it is of some interest to examine the predic- 
tions of singular perturbation theory in this context. 

To do this, we consider a simple model with a scalar 
NCOP - the Fisher Equation - and add a convective 
nonlinearity, whose strength can be tuned with a dimen- 
sionless parameter /x. We call this new equation FEC 
(for Fisher Equation with Convective nonlinearity; cf. 



eq.(l)). The Fisher equation has been extensively stud- 
iedin literature [7-9], originally in the context of pop- 
ulation dynamics, and the KYG analysis yielded good 
results from its solution [7](c/. end of Introduction). For 
FEC with /X < 1, the KYG SPM gives results very sim- 
ilar to those found for the Fisher equation, but breaks 
down at early times when /x ~ 0(1) or greater. This 
is because when /x is large enough, mechanical transport 
(which is a non-linear process) dominates over diffusive 
processes (which are linear). The KYG analysis can be 
expected to give usable results only when the front veloc- 
ity is selected by the linear dynamics. This effect of the 
convective nonlinearity is expected to be general to other 
systems where bias field or convection is important. 

Before we proceed to describe our results, some gen- 
eral remarks about the validity and utility of the KYG 
technique are in order. This technique is not meant to 
be a means of obtaining a solution to the initial-value 
problem for a reaction-diffusion equation. Rather, the 
singular perturbation method should be interpreted as 
a way of obtaining an analytic form, which may repli- 
cate the important features of the true solution, thereby 
enabling one to obtain statistically important quantities. 

KYG [1] first applied the singular perturbation method 
to the Time-Dependent Ginzburg-Landau (TDGL) equa- 
tion with a scalar order parameter. In this case, the ana- 
lytic solution obtained is in disagreement on one impor- 
tant point with the real (obtained numerically) solution 
of the TDGL equation [7], viz. it has infinitesimally thin 
walls at late times, whereas the real solution always has 
walls of nonzero thickness. However, the analytic solu- 
tion does reasonably reproduce the defect (interface) dis- 
tribution of the real solution, starting from random initial 
conditions. Also the time-dependent structure factor cal- 
culated from the analytic solution of KYG [1] is identical 
to the better-known result of Ohta et al. [10], which is 
derived using interface dynamics. Thus, for the scalar 
TDGL equation, the singular perturbation approach is 
able to predict the "correct" asymptotic structure factor 
and the growth law for the characteristic domain size ( 
~ t^/^, where t is time). The next application of the 
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SPM was to the d-dimensional Fisher equation [7,11]. 
As in the TDGL case, the interfacial profile of the KYG 
analytic solution (which is a travelling wave front) is ar- 
tificially sharp compared to that for the real solution. 
However, the (travelling wave) analytic solution has the 
correct asymptotic velocity, even though the approach to 
this asymptotic velocity is incorrect in one dimension, 
where the exact result is known [9]. In the case where 
long-range interactions are also present, these statistics 
and defect dynamics are actually incorrectly given by the 
KYG method [12]. We find similar limitation for our non- 
conserved order-parameter system described by the FEC 
equation, in contrast to other non-conserved systems. 

In Section II, we extend our earlier applications of the 
KYG method [2,5,7,11], to the FEC, and follow the for- 
malism described in detail by Puri [7]. Because the KYG 
analysis can be applied to FEC with only small differ- 
ences as compared to that done for the Fisher equation, 
only an outline of the method is given here, emphasising 
these differences, and the reader is refered to [7] and [1] 
for further details. In Section III, we compare the ana- 
lytical solution with numerical results and summarize. 

II. MODEL: KYG METHOD ON FEC 

The FEC in one space dimension is 

du du d^u 2 

where the order parameter field u, the position x and the 
time t have been rescaled to dimensionless units. The 
dimensionless quantity /x is the ratio of convective to dif- 
fusive strengths in the equation of motion for u and can 
not be scaled out. This equation is equivalent to Burg- 
ers equation [13] in which one would add a linear source 
term u (providing the instability) and a quadratic sink 
term (providing the damping). Burgers equation is 
the simplest PDE to combine a nonlinear convective term 
with linear diffusion and is analytically solvable. It often 
appears as a limiting case in more complicated hydro- 
dynamics problems, usually in the context of turbulence, 
and its main feature, due to the convective term, is that it 
yields shocks, also observed in the FEC. Note that /x = 
is the Fisher equation itself. 

The essential idea of the KYG method is to gener- 
ate an infinite order perturbative expansion around the 
linear solution of the PDE, approximate the n*'' order 
term in a suitable manner and resum the resulting in- 
finite series to get an approximate, yet analytic, closed 
form result. There is an alternate derivation (also by 
KYG) using functional analysis that can be used for sys- 
tems in which the non-linearities do not contain gradient 
terms. Since this is not the case here, we stick to the 
"standard" method. 

The solution, in Fourier space, of the linearized equa- 
tion (1) is ul{t) = eT'=*'Ujt(0) with = I - ■ Thus for 



k < 1, ■u°(i) grows exponentially in time with rate jk- 
In what follows, the tilde denotes Fourier space and the 
superscript denotes the solution to the linearized PDE. 

The nonlinear equation (1) can be rewritten in an in- 
tegral form in terms of the linear solution as 




g{h)uk,{t')uk,{t')dhdk2dt' . (2) 



where g[k) = 1 — ifj,k and 6 represents a Dirac delta 
function. The effect of convection is totally included in 
g{k), which would be identically 1 for the Fisher equa- 
tion. In applying the KYG technique, this modification 
creates no difficulty. Note however that if the nonlinear 
damping term had been a instead of a , g[k) would 
actually be a functional of u, so that the KYG analysis 
would be impossible to carry out. 

We consider the integral in (2) as a singular per- 
turbation: we introduce an expansion coefficient A in 
front of it and expand Uk{t) in a power series in A, 
Uj.{t) = "^n-Q An{k,t)X^ . By equating the coefficients 
of equal powers of A, expressions can be found for the 
An{k,t) in terms of ■u°(i), and the result put in a dia- 
grammatic form. Each An grows exponentially faster in 
time than An-i- Higher order An rapidly get hopelessly 
complicated so that each of the n\ diagrams making up 
a given An{k,t) is approximated by the "comb" diagram 
(denoted C^^'^ [k,t)) of the same order. The absolute er- 
ror made in this approximation (which essentially gives 
equal weight to two- and multiple-A;-mode time correla- 
tions) is uncontrolled and constitutes the most drastic 
approximation in the method. 

At this point we have 




which is identical to the form obtained for the Fisher 
equation except for the factors 3(^:2^-1) in the integral. 
Using a Laplace transform, (3) can be separated into n+1 
terms, each growing exponentially in time. In order to 
simplify and to keep the expression tractable, a second 
approximation is made by KYG whereby only the fastest 
growing term of the sum is kept, which induces a relative 
error exponentially decreasing in time. This approxima- 
tion is not affected by the form of g[k) as long as g[k) is 
independent of Uk{t) (which is the case for FEC). 
The integrand, now, is of the form 

/({A=J)ff(A=2.-i)e''«'='»ii°^,_^(0)ii°^^(0), 
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where / and h are functions whose form are not impor- 
tant for now, and where {ki} denotes the set of ki's 
appearing in the integral. At late times one expects 
the maxima of h to dominate the integrand, so one 
can apply the Laplace integral method to simplify it. 
One can show that the exponential is maximized when 
k2i-i = k/[n+ l),Vi. There are several different ways 
of approximating the integrand, e.g. one can choose to 
evaluate only /({^i}) at the maximum, or both /({^i}) 
and g[k2i-i), but only one yields a final result which is 
consistent with the initial PDE: one must evaluate both 
/({^i}) g{k2i-i) at the maximum of the exponential. 
Applying these various steps to (3) yields: 
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This means CP'^ = o(-^ 



■yFE 



[9'(,;r+T-'J <"n"' ^^"^ super- 
scripts FEC and FE refer to the Fisher Equation with 
Convective nonlinearity and the Fisher Equation results 
for C„, respectively. Even with the much simplified 
comb diagram (4), the power series for Uk[t) is still not 
summable analytically. To do this and also be able to 
Fourier transform back to real space, one usually approx- 
imates the coefficient of the multiple integrals in (4), 
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to first order in k, since the major contribution to the 
Fourier transform comes from k near (for NCOP sys- 
tems only). For /x = this is a very good approximation. 
But for /X 0, g[k) is complex, so that a phase error is 
also introduced. Because our singular perturbation ex- 
pansion contains an infinite number of terms which must 
partly cancel each other to give a convergent sum, this 
phase error, although small, could introduce cancellation 
(interference) effects whose consequences can be judged 
only a posteriori. 

With these approximations, one can analytically resum 
the power series in A and invert the Fourier transform to 
get [15] 



FEC/ j\ 
U [x, t) 
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Note that we recover the KYG solution to the Fisher 
equation when /x = 0. The fj,^ term induces the same 
type of asymmetry in (6) as that given to the PDE (1) by 
the convective term. However one of the consequences of 
fj,-^ is the development, at large times, of strong peaks 
in u[x,t) at the interfaces, indicating a breakdown of 
the singular perturbation result. Furthermore, one can 
see that this occurs for any finite order expansion of (5) 
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FIG. 1. Numerical solution of FEC at i = 20 and sev- 
eral values of \i [cf. figure legend), from a localised ini- 
tial-condition in a; = 0. The /i = curve is symmetric with 
respect to the origin. Others get distorted by the convection 
term [/i ^ 0). 



in k. Given that higher order diagrams [n — > co) will 
dominate the power series for U}.(^) at large times, and 

that the k values dominating the dynamics are near 0, a 

r 1 ^ 

remedy would be to approximate g'(;;;^) with a sim- 

r 1 ^ 

pier non-polynomial function in k, close to g'(;;;^) 

for A; — > and n — > co. The one possibility which will 
give tractable summation and a fourier-transformable di- 
agram is exp(— i/xA;). This approximation, which is very 
good for /X <C 1 (c/. next section) exhibits no peaks at 
the interfaces even as i — > co. Then from (4) one now 
gets, instead of (6), 



PEC, 



E, i) ~ exp(/x 



dx ' 



u°{x,t) 



(7) 



Although the exp(/x-^) operator does not create ar- 
tificial peaks at the interfaces for any /x and t, it is in 
effect a translation operator. Hence if the KYG solution 
to the Fisher equation (/x = 0) is denoted by u^^[x,t) = 
u°{x,t)/{l + u°{x,t)), we have u^^^{x,t) = ■u^^(a;-/x, i), 
i.e., the solution to the FEC we get from the singular per- 
turbation expansion is in fact a translated solution to the 
Fisher equation. Since the true effect of the convective 
term can not be equated to a translation, it is immedi- 
ately apparent that important characteristics of convec- 
tion have been lost in this singular perturbation analysis. 
This will be crucial when /x is sufficiently large. We now 
briefiy discuss and compare these analytic results with 
those of numerical integrations. 



III. DISCUSSION 

In Fig. 1 we show the exact numerical integration of 
the FEC equation, at relatively small values of /x and 



3 



intermediate times [t = 20), for a seed initial condition 
centered at a; = (i.e. u[x, 0) = 6[x)). The /x = curve 
is a numerical solution to the Fisher equation. The stan- 
dard Euler integration method was used for all numerical 
integrations, with time and space meshes suitably chosen 
to insure stability and precision. After a transient time 
approximately equal to 10, the interface profiles change 
little and the interface velocity reaches a value within 
10% of its asymptotic value. Hence the curves for i > 20 
are almost identical to those of Fig. 1 but with larger 
bulk region (where u ~ 1). As for the shock waves found 
in the solution of the Burgers equation [16], the effect of 
the convective term on the spatially symmetric solution 
of the Fisher equation is to steepen the right interface 
and broaden the left one, as /x increases. One way of see- 
ing analytically how the convective nonlinearity operates 
is to look at the PDE -^u = which has, in implicit 

form, u[x,t) = f[x — u[x,t)t) as a solution. This is a 
"wave" with points at height u travelling at speed u, so 
that higher points travel faster than lower ones. 

In Fig. 2 we compare the results of a numerical inte- 
gration for FEC with the analytical results of equations 
(6) and (7), for /x = 0.1 and t = 20. For such a small /x 
the asymmetry is difficult to discern by eye for all three 
curves, but can be checked numerically. Also, because of 
the smallness of /x, it is hard to distinguish between the 
solutions given by (6) and (7) (as expected from the ap- 
proximations), so that one could use either (6) or (7) as 
approximate solutions to FEC, keeping in mind the re- 
marks made at the end of last Section. The offset between 
the interface of either of the approximate solutions and 
that of the exact numerical one occurs because at early 
times the speed with which the fronts in the exact nu- 
merical solution travel is significantly smaller than that 
of the approximate solutions: the velocity of the fronts 
given by (6) and (7) are both equal to the asymptotic 
velocity of the exact solution, namely 2, as for the Fisher 
equation [7]. 

The interfaces given by our approximate solu- 
tions can be softened further by approximating 

nr=i (l + {i+Tr''^) in (4) by exp{-P/2) instead of 1, 
which is an extremely good approximation (less than 1% 
error) for all A; < 1 and n > 4 (this softening is analogous 
to one done by Oono and Puri [17] and Puri [7]). How- 
ever one finds that the change this induces in the profiles 
is very slight (not shown for that reason) and hence only 
useful at early times, so that the hardness of the inter- 
faces of the approximate solutions is more deeply buried 
in the comb diagram approximation. 

For /X ~ 1, the asymmetry between the two sides of 
the solution is much more obvious. But also (as shown 
in Fig. 3) eq. (7) gives a better solution than eq. (6). The 
steepness of the right interface of the eq. (6) solution at 
i = 40 indicates the beginning of the breakdown of the 
SPM result, since going to slightly larger times brings 
about a clear peak in the left interface and a dip in the 
right one. This does not happen for eq. (7). From the 
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FIG. 2. Order parameter profile for the numerical solution 
(solid curve), for eq. (6) (dotted curve), and for eq. (7) (dashed 
curve), at i = 20 and /i = 0.1. 
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FIG. 3. Left and right interface profiles of the numerical 
solution (solid curve), of eq. (6) (dotted curve) and of eq. (7) 
(dashed curve), at i = 40 and fi = 1. 

point of view of defect dynamics, the two solutions give 
the right velocity of 2 (obtainable from linear dynamics). 

For large /x (order of 10), both approximate solutions 
break down. The asymmetry of the analytical solution is 
clear but it matches in no way that of the real numerical 
solution. Indeed, the velocity of the right interface is 
wrongly given by both analytic solutions, as it is much 
greater than 2 (c/. Fig. 4) at all times. Peaks develop 
at the interfaces given by eq. (6) even at intermediate 
times like i ~ 20 (not shown). Also, the left interface of 
the exact numerical result is extremely broad for /x on 
that order, something that the SPM seems incapable of 
describing. This can hardly be taken into account by the 
approximations used and, as can be seen, both eqs. (6) 
and (7) break down for the relative convective strength 
of the order of 10 (it actually breaks down for /x ~ 5 at 
t ~ 40), not only in terms of order-parameter profile, but 
of defect dynamics. 
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or convection is present. 
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FIG. 4. Numerical solution (solid curve) and eq. (7) 
(dashed curve) at i = 20 and fi = 20. Again, from initial 
condition in a; = 0. 

We can therefore say that KYG SPM gives results very 
close to the real solution in terms of profile and velocity 
up to /X ~ 0(1). But as the shock wave nature of the 
convective nonlinearity manifests itself more and more 
strongly (as /x increases), the KYG SPM does not man- 
age to capture the essential features of convection, both 
in profiles and defect dynamics. Several analytical ap- 
proaches are being investigated by us to make the results 
more quantitative and insightful. Apparently there is 
more to the failure of KYG at large /x than just a limita- 
tion of the SPM: any type of analysis performed around 
the linear solution, as done here, cannot be expected to 
work if the (attracting fixed point) solution becomes hy- 
perbolic at some large enough value of /x, while an al- 
ternate solution (corresponding to a different asymptotic 
profile and velocity) becomes stable. 

IV. CONCLUSION 

The utility of singular perturbation methods lies in the 
calculation of statistical quantities (e.g., time-dependent 
structure factors, domain growth laws), which are deter- 
mined by the qualitative features of the solution. How- 
ever we have shown that the singular perturbation ap- 
proach will not give a reasonable solution to the initial- 
value problem for a reaction-diffusion equation where 
convection is strong enough (a precise criterion will be 
given in a future article) even in terms of such statisti- 
cal quantities since the defect dynamics are incorrectly 
predicted. It can be shown [18] that this incorrect defect 
dynamics is due to the convection term dominating the 
diffusive transport of perturbations, and is not due only 
to some intrinsic limitations of the singular perturbation 
approach itself. Hence approaches which treat the full 
non-linear PDE must be sought to study the dynamics 
of reaction-diffusion systems where an external bias field 
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